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Abstract 

We examine a general multi-factor model for commodity spot prices and futures valuation. We extend 
the multi-factor long-short model in [T] and [2] in two important aspects: firstly we allow for both the 
long and short term dynamic factors to be mean reverting incorporating stochastic volatility factors and 
secondly we develop an additive structural seasonality model. Then a Milstein discretized non-linear 
stochastic volatility state space representation for the model is developed which allows for futures and 
options contracts in the observation equation. We then develop numerical methodology based on an 
advanced Sequential Monte Carlo algorithm utilising Particle Markov chain Monte Carlo to perform 
calibration of the model jointly with the filtering of the latent processes for the long-short dynamics 
and volatility factors. In this regard we explore and develop a novel methodology based on an adaptive 
Rao-Blackwellised version of the Particle Markov chain Monte Carlo methodology. In doing this we 
deal accurately with the non-linearities in the state-space model which are therefore introduced into the 
filtering framework. We perform analysis on synthetic and real data for oil commodities. 

Key words: Multi-Factor, Commodity, Spot Price, Stochastic Volatility, Milstein, Adaptive Markov 
chain Monte Carlo, Particle filter, Rao-Blackwellization. 



1 Commodities economic theory 

The basis of commodity modelling has a long history in the economics and econometrics literature, indeed 
instrumental concepts were introduced in [3j which developed normal backwardation to attempt to model 
inter-temporal futures prices, in order to explain differences between spot and futures prices. Storage theory 



was developed in [1] to account for storage costs, which basically states that inter-temporal prices are the 
jointly determined price of storage. That is, futures prices will be approximately the spot price plus positive 
storage costs. One problem with this explanation was that the interpretation of backwardation, which 
occurs when the observed Futures price is less than the Spot price, was difficult to explain. To address 
this problem with interpretation, the work of [5j developed the concept of convenience yield, which is the 
implicit revenue associated with stock holding. This can be considered as analogous to coupons linked with 
bonds or dividends payed on stocks. Put another way the convenience yield can be interpreted as the flow 
of services accruing to the holder of the spot commodity but not to the owner of a futures contract. 

The conclusions of storage theory provided the following fundamental elements that models of commodi- 
ties should aim to address when considering the economic relationship between spot and futures prices 
implied by storage cost. This relationship identifies at least three key factors influencing futures prices, 
given by the Spot price, the convenience yield and the interest rate. Additionally, convenience yield and 
spot price are positively correlated and are both inverse functions of stock level. One also needs to consider 
that there is an asymmetric relationship for basis, that is the variation between spot price of a deliverable 
commodity and the relative price of the futures contract closest to maturity. In particular when the market 
is in contango (futures > spot), the basis level is limited to the storage costs. One final important economic 
relationship that commodity models should aim to address is the Samuelson Effect, [6], which states that 
movement in prices of prompt (soon to mature) contracts are large and erratic whilst prices of long-term 
contracts are not erratic. This implies that variance of futures prices and correlation between nearest futures 
price and subsequent contract prices decrease as a function of time. 

1.1 Outline and Contributions 

In this paper we begin with a review of multi-factor commodity futures models. Then in Section 3 we 
propose and develop a general multi-factor model for commodity spot prices and futures valuation. We 
extend the multi-factor long-short model in [Ij and [2] in two important aspects: firstly we allow for both 
the long and short term dynamic factors to be mean reverting incorporating stochastic volatility factors 
and secondly we develop an additive structural seasonality model. In Section 4 a Milstein discretized non- 
linear stochastic volatility state space representation for the model is developed which allows for futures and 
options contracts in the observation equation. In Section 5 we develop an Bayesian filtering and calibration 
framework in which we develop numerical methodology based on an advanced Sequential Monte Carlo 
algorithm utilising Particle Markov chain Monte Carlo to perform calibration of the model jointly with 
the filtering of the latent processes for the long-short dynamics and volatility factors. Section 6 explains 
this novel methodology in detail. In this regard we explore and develop a novel methodology based on an 
adaptive Rao-Blackwellised version of the Particle Markov chain Monte Carlo methodology. In doing this 
we deal accurately with the non-linearities in the state-space model which are therefore introduced into the 
filtering framework. In Section 7 we perform analysis on synthetic and real data for oil commodities. 



2 Review of multi-factor dynamic commodity futures models 



The typical approach to formulating a model for commodities is to begin by specifying a real world latent 
stochastic differential equation (s.d.e.) based model for the factors to be considered in the commodities 
model. Then under assumptions on the market prices of risk and the behavior of such markets, a change 
of measure is made to adjust the real world process by the price of risk in order to obtain the risk neutral 
measure. Under this risk neutral s.d.e. model, the futures price are obtained as functions of the latent 
process factors as expected discounted future payoffs. 

In [7j a reduced form model for commodity prices was introduced. This was a one-factor model in which 
the commodity price follows a Geometric Brownian motion (GBM) and the convenience yield is treated as 
a dividend yield. This model was clearly not sufficient as it failed in two important aspects, firstly it does 
not account for the property of mean reversion of spot prices and secondly it neglects inventory dependence 
properties of the convenience yield. This basic model was extended by the work in [8] to a two factor model 
and used to study oil futures. In the two-factor model the spot price was modelled by Geometric Brownian 
Motion and the Convenience Yield was modelled by a mean reverting process. 

As noted in [T], this two factor model can be made more parsimonious via a change of variable. A new 
state variable is defined by the demeaned convenience yield minus the long-term convenience yield and at 
time t this represents the short-term deviation in log prices or the deviation of spot price returns from its 
long-term value. Next, define a second factor as the long-term price return or price appreciation which is 
given by the long-term total return minus the long term convenience yield. This reformulation in addition to 
being parsimonious is also easier to interpret from a physical perspective since now the returns are defined 
in terms of the long-term price appreciation. 

In the three- factor model, ([9]; [10]; [Tl]) a third unobserved stochastic variable is introduced, with the 
short term variation (follows OU process) , the equilibrium price (follows GBM process) and the instantaneous 
growth in equilibrium prices (follows a mean reverting process) for the real process. 

Note, it should be pointed out that an additional assumption on these models is that there are no 
arbitrage opportunities. Strictly speaking this will not be true since, as pointed out in [12], these models do 
not guarantee that the convenience yield is strictly positive. This is required for these models to be arbitrage 
free, where contemporaneous spot prices are greater than discounted futures prices net of carrying costs. 
Additionally, as pointed out in the introduction, a key element of storage theory is that the variance and 
correlation between spot and futures prices should be related to the level of the price or convenience yield. 
These models make the simplistic assumption of constant correlation and variance. Finally, in addition to 
these weaknesses, the spot price volatilities obtained from these models are strongly hetroskedastic, see [12j 
for a review of literature on studies demonstrating these issues. 

More recent models include those developed in [12] which removes the potential for cash and carry 
arbitrage, by forcing the convenience yield to be positive at all times. Additionally, these new models admit 
time varying spot and convenience yield volatilities. This model modifies the models discussed whilst still 



admitting an analytic expression to be obtained for the futures contract price as a function of the unobserved 
stochastic factors. In particular the two factor model they develop has a spot price following a GBM in 
which the convenience yield is treated as an exogenous dividend yield and the volatility is proportional to the 
square root of the instantaneous convenience yield level. The second factor, the convenience yield, follows 
a Cox-Ingersoll-Ross (CIR) model, forcing convenience yield to be strictly non-negative, with volatility 
proportional to the square root of the instantaneous convenience yield level. Under the assumption that 
the interest rate is constant, an analytic solution is obtained for the futures price as a function of these 
unobserved stochastic factors (real process). 

Another recent model is from [2j which developed a very general multi-factor model incorporating several 
extensions to the above models including a Levy driven process with jumps. However, again in this model 
the questions relating to joint parameter estimation and filtering were not addressed, nor was seasonality 
incorporated into the model. The stochastic model developed was given by (risk neutral setting): 



Where St is the spot price at time t, 6t is the instantaneous convenience yield at time t, the long term 
total return (price appreciation + convenience yield), k^, Ky and ks are the mean reverting coefficients, Or 
and 9v are the long term interest rate and volatility, as,(Tr the volatilities of spot price and convenience 
yield dynamics. 

We modify this model structure for this paper and derive the resulting futures price as a function of 
the latent commodity factors. To achieve this we remain in the exponentially affine family of models. The 
derivation of the futures price for our model allows us to use it as the observation equations for these models 
that we consider in this paper which is derived as a exponentially affine linear function of the state variables. 
The observation equations linking these unobserved factors to the futures price were obtained under the risk 
neutral measure. 

3 Proposed Multi-Factor Stochastic long-short price dynamics model 
incorporating seasonality. 

In this paper we will extend the long-short stochastic Model to allow both the long and short stochastic 
factors to be mean reverting processes and we incorporate an additive structural model to account for 
seasonality and a third factor for stochastic volatility to allow for modelling of implied volatility smiles. 



dSt = in -6t- Xfij) Stdt + asStdWi + 





drt = {Or — i^r-ft) dt + ar^/rtdWr^, d6t 



{9s - n&5t) dt + asdWs 




3.1 Latent stochastic factors and seasonality. 

Here we present the model considered in the remainder of this paper. 



Real World Process 


Risk Neutral Process 


X{t)=ln{St) = Xt + ^t + 0t + f{t) 


X* (t) 


= X*t+Q + 0* + f{t) 


dxt = -Pxtdt + a^dZ^ 


dx*t = 


{-/3*X*t - Ao) dt + a^dZ*^ 


d^t = (/U? - i^s,it) dt + a^dZ^ 


dCt = 


(^1 - i^l^t) dt + a^dZl 


dOt = ^/VtdZe 


det = 


-XiVt*dt + ^/WdZ; 


dVt = (/iy — KyVt) dt + avVVtdZv 


dVt* = 


itJ,*y - K*yVn dt + ay^dZ-y, 



with correlations structures given by: 

E \dZ^dZQ\ = E [dZ^dZg] = E [dZ^dZy] = E [dZ^dZy] = 0; E [dZ^dZ^] = p^^,; E [dZydZg] = pyg. 

In addition we define f{t) as the deterministic unknown seasonal component and the risk premia in the risk 
neutral model are given by: 

Ax {Xt,it, Bu Vt) = \l + Kxu h ixt, Ct, Ot, Vt) = X*2 + A^6 

Ae {xt,Ct,Ot,Vt) = \l + Xl^t; Xv {xt,^t,et,Vt) = Xl + X;Vt 

Ao = AqC^x; Ai = Xla^; X2 = X2(J^; A3 = X^a^; A4 = A4; A5 = A5 = 

Ae = Agcry; A7 = Ayciy;/?* = /3 + Ai; = /i^ - A2 

k| = Kg + A3; py = i^v - Ae; Ky = ny + At- 

Note the real processes are simply obtained by setting the risk premia (A) parameters to zero. We can 
also ensure this model remains arbitrage free, according to discussions in [13], by considering constraints 
on the parameter py such that 2py > 1 under both the real world and risk neutral measures, to ensure 
boundary non-attainment discussed in [13j. This can be enforced in the prior specification we develop in 
Section 5. The parametric seasonal model we considered in this paper is given by: 

12 

fit) = ^ujknkt, (3.1) 

k=2 

where ^Ikt = 1 if t is the date of the k-th calendar month and ^Ikt = otherwise. Note, without loss of 
generality one can set /(O) = as this constant is absorbed by and all variation is with respect to January. 

3.2 Commodity futures price. 

In this section we obtain analytic expressions for the futures price under the proposed multi-factor long- 
short dynamic model. The futures price we consider here will be developed as a function of the underlying 
stochastic factors. We denote the Futures contract price at time t with maturity of the contract at T by 
Ft^T with the proof of this results provided in Appendix 1. 

ProDOsition 1: The futures contract vrice F+t at time t for a contract maturina at time T has a closed 



from solution, with respect to the latent state variables that comprise the log spot price, Xtj^tT^t) when priced 
under the risk neutral measure, given by 

Ft,T = exp (/(T) + Bo{t) + Bi{T)Ct + B2{T)xt + B^{T)et) , 
with coefficients given by, 

2 2 

I exp(-4r) - ^ exp(-rr) + expC-^r - ^r) 

Bi{t) = exp(-4r);52(T) = exp(-/3V); 53(t) = 1. 

The proof of this result is provided in detail in Appendix 1. 

Remark 1: Note that the parameters Bi{T) are under the risk neutral framework and the latent factors in 
the futures price are real. 

Having derived the futures price for the stochastic model we propose, we now specify a state space model 
formulation and a Bayesian model. We will then develop a novel sampling methodology to estimate jointly 
the model parameters and the latent factors. 

4 State Space Model via Strong Milstein Sheme: long-short stochastic 
volatility model. 

In specifying the state space models, we shall consider a panel of N different futures contracts with maturities 
Ti, . . . , Tat and futures prices at time t denoted by -Ft,Ti, • • • , F't^TN- In addition, we note that the filtering 
framework is set up such that the latent process is the real process with real parameters used and the 
observation process is the real process with risk neutral parameters used. 

4.1 State Equations 

In this section we convert our s.d.e. models for the latent factors into a discretized version to obtain a 
state space model formulation. To achieve this we consider strong Taylor schemes, where an approximation 
process Y converges in a strong sense with order 7 G (0, 00] with a continuous process X if there exists a 
finite constant K and a positive constant 6 such that, 

E{\Xt-Yn\) ^ K6^. 

We employ a strong stochastic Milstein Scheme. For the long-short factors the strong stochastic Taylor 
scheme of order 1, otherwise known as a Milstein scheme, corresponds to the Euler scheme. 

For the volatility related stochastic processes the Milstein scheme is employed for these two processes 
since firstlv thev are both correlated, so the same discretization scheme should be emoloved for each of 



them. Secondly the process for dVt certainly does not have a Gaussian transition density, as such a Euler 
type scheme with additive Gaussian increments is too simplistic an approximation, even at daily intervals. 
In the Milstein scheme we truncate the mixed integrals at a p-th order as shown in Appendix 2. 

Proposition 2: The bivariate s.d.e. 's for xt and S^t with 2- dimensional Wiener process discretized under 
a strong Euler scheme leads to the following state space equations. 

Xt = il- /3At)xt_i + ^xV^i^x 

6 = + (1 - K^At) ^t^i + a^\/Mn^ 

where n^,n^ are standard normal random variables for the innovation noise. The bivariate s.d.e. 's for Ot 
and Vt with 2-dimensional Wiener process discretized under a strong Milstein scheme leads to the following 
state space equations. 

Ot ~ Ot-i + y^Vt-iAtng^t + ^avpve {'^tnj ^ - At) + ^crvyjl - Pve'^^2,i)At 

Vt w Vt-i + {pv - KyVt^i) At + avpvey^Vt-iAtng^t + crv\JVt-i (l - Pyg) Atnv,t 

+ \(^vPve (Atn^,t - At) + ^(^l^pve^Jl - ple'^^i,2)At + ^dym-^l - Jg ^^^^ + (l - Pyg) (Atn^^ - At) 

where ng^t and ny^t are i.i.d. standard normal random variables and we utilise a p-th order truncation with 
terms for J^. .^^^ defined by 

'^Uuh)At = (\ChCj2 +VPp(.Ph,pCj2 - Pj2,pCh)) (^ii.^ (^^2 + l'j2,r) - i'j2,r [^Cji + I'jur)) , 

\ J 

where Cj^ Pj,p:'^j,r andipj^r are all independent M {0; 1) Gaussian random variables with, 

1 1 -A 1 1 

r=l ^ 

In the specification of our models we take At = 1 day. 
The proof of these results is presented in Appendix 2. 

Remark 2: Under this discretization scheme it will not be possible in general to find the transition 
density for these two factors 6^ and Vt, given 6f^i and Vt-i. We note that the version of the particle filter 
we develop only requires that we can generate from these state equations, which is trivial in the sense that 
it only requires generating Gaussian random variables and applying a transform. This is the key reason for 
our use later on of the SIR particle filter approach in our estimation methodology. 

4.2 Observation Equations 

Clearly, the observation equations for the futures contracts are linear and will be given in Equations 15.11 
lnFi,Ti = /(Ti) + Boin) + Bi{n)Ct + B2{n)xt + Bs{n)9t + ^ 



In Ft^Tj, = f{TN) + Bo{tn) + Bi{TN)^t + B2{TN)Xt + Bs{TN)et + el 



where e\ ' , . . . ,el ' are all zero mean normal random variables for the observation noise on the futures 
contract prices, with covariance structures at each time t denoted by which is {N x N) dimensional. 
Clearly in this model the state equations, presented in Proposition 2, are non-linear and non-Gaussian, 

hence the use of a Kalman-Filter is no-longer optimal. It is also clear that a simple linearization of the state 
equations will not suffice, hence this also precludes the use of the approximate Extended Kalman Filter 
formulation. Instead we will develop a particle filter and Markov chain Monte Carlo solution based around 
a novel version of the Particle MCMC algorithm. 



5 Bayesian Model for Commodities 

In this section we develop a Bayesian model for the parameters and latent factors that must be estimated 
over time. In introducing a Bayesian formulation for this problem we note that we model the real process 
parameters as random variables and the risk premia used to convert to the risk neutral pricing framework as 
unknown random variables. The Bayesian formulation proceeds by first defining the posterior distribution 
one is interested in considering: 

TT (*, Xl:U ^1:U 01:t, '^^1:* |-Pl:t,Ti , • • ■ , Pi.t^rJ (5.1) 
T 

= /x*(xo, ^0, ^0, Vo) n /# ^ (XtlXt-i) /^^ mt-i) {0t\et-i, V^-i) fS'^ {Vt\Vt-i) (5.2) 

T N 

^T[T[9'^{Ft,T^Xuit,9t)pm (5.3) 

t=l n=l 

where we define quantities as follows: 

• $ = ^jjjV) "^iiii)) where and ^rn Sive the real and risk neutral parameter vectors respec- 
tively, e.g. = (^, /Lt^, K^, /Lty, Ky, (T;^, (T^, (Ty, Eg) the random vector of unknown real-world model 
parameters. 

• A priori there is no reason to specify the priors for the real-world and risk-neutral parameters differ- 
ently and therefore we give them identical prior specifications. With the real-world prior p(^r) for 
the random parameter vector ^r given by: 

With respect to the choices of hyper-prior parameters, we consider an objective uninformative prior 
specification, see Section 7. Furthermore, we assume that = diag (i^g^, . . . )0"ejv)' though it is trivial 
to consider alternative specifications. 

• ;^#(xO) ^0) ^0) ^o) is a multivariate Gaussian with known mean u and covariance Tjh which specifies the 
distribution of the initial state of the factors eiven the model Darameters. 



• {Xt\Xt-i) is given by the state space model state equation to be a Gaussian with mean (1 — 
/3)Xt-iAt and standard deviation u^-^/St. 

• (CtlCt-i) is given by the state space model state equation to be a Gaussian with mean ^u^At + (1 — 
K^At)xt-i and standard deviation a^^^J /\t. 

• {Ot\dt-i) and (V(|T4_i) are given by the state space model state equations to be the distribu- 
tion formed from the addition of the random variables ^/Vt-ine t + + e^'^ and av\/Vt~inv,t + 
respectively. The convolution integrals involved with finding these distributions do not admit a para- 
metric form. We that this does not affect our estimation approach since we only require that we can 
sample efficiently these processes. 

• g^{Ft^Tn \xt-,it-, dt) is given by the state space model observation equation to be a Gaussian with mean 
/(Ti) + Bo{t) + Bi{t% + B2{T)xt + Bz{T)et and standard deviation 1. 

6 Parameter estimation and filtering via Adaptive MCMC within Rao- 
Blackwellised Particle MCMC 

In this section we detail the simulation methodology we develop for obtaining samples from the posterior 
distribution given in Equation (15. 3p . This simulation methodology is based on a recent sampler specifically 
developed for use in state space models in which one wishes to estimate jointly the model parameters and 
the latent states, see [Hj. The methodology is known as Particle Markov chain Monte Carlo (PMCMC) 
and represents a state of the art sampling framework for such problems. The approach embeds a particle 
filter estimate of an optimal proposal distribution into a MCMC algorithm, allowing one to obtain samples 
from the target posterior distribution for the parameters and latent states. This embedding is done in a 
non-trivial fashion, see [Tl] for details. 

The Particle MCMC proposal distribution is split into two components. The first involves a proposal 
kernel which is constructed via an adaptive Metropolis scheme and this will be used to sample the static 
parameters The second component of the proposal kernel involves the sampling of a trajectory for the 
latent factors, Xi:t) This proposal kernel is constructed via a Rao-Blackwellized Sequential 

Monte Carlo algorithm in which the Rao-Blackwellizing variance reduction is achieved via optimal Kalman- 
Filtering for the path space associated to Xi-.tiii-.t- 

The introduction of an adaptive MCMC proposal kernel into the Particle MCMC setting allows the 
Markov chain proposal distribution to adaptively learn the regions in which the marginal posterior distri- 
bution for the static model parameters has most mass. As such the probability of acceptance under such 
an online adaptive proposal will be significantly improved over time. Then for the latent factor path spaces 
we develop a minimum variance approach involving Rao-Blackwellisation via a Kalman Filter, followed by 
a Sampling Importance Resampling (SIR) version of the particle filter. Note - it is important to work with 



an SIR filter for the path space parameters associated with 9i.t,Vi-t. The reason for this is that under the 
Milstein discretization scheme we develop, we can efficiently sample exactly from the prior parameters but 
can not evaluate the prior density point-wise. Hence, under this variance reduction sampling framework, we 
optimally sample the latent state path space trajectories for xi-.t and via a Kalman Filter formulation 
and we use a SIR version of the particle filter for the state path trajectories for 6i-t, Vi:t. 

6.1 Particle Markov chain Monte Carlo (PMCMC) Background 

Particle MCMC is a new methodology introduced recently in [14J . The methodology utilizes particle filtering 
methodology to approximate the optimal proposal distribution in an MCMC algorithm in which one is 
interested in high-dimensional block updates. This naturally, lends it self well to the setting of state space 
modelling. In the setting considered in this paper the parameter space of static parameters in the model is 
high dimensional, and in addition the path space latent factors given an additional 4t parameters. Where, 
in this case t can range from hundreds to thousands of days, depending on the commodities studied and 
sampling frequency of the observations. Hence, this problem is ideally suited for the PMCMC methodology. 

The simplest solution to this problem would be to sample directly form the univariate full conditional 
distributions. However, even in the case in which a single component Gibbs sampler is achievable via 
inversion sampling for all parameters, this sampling scheme will typically results in very slow mixing of 
the Markov chain around the support of the posterior. This is especially problematic in high dimensional 
target posterior distributions, since it requires excessively long Markov chain runs to achieve samples from 
stationary regime. Typically, it can also lead to very high autocorrelations in the Markov chain states. 

It is well known that to avoid this slow mixing Markov chain setting, one must sample from larger blocks 
of parameters. However, the design of an optimal proposal distribution for large blocks of parameters is very 
complicated. The Particle MCMC methodology allows one to approximate the optimal proposal distribution 
for a large number of parameters via a Sequential Monte Carlo approximated proposal distribution. 

In the state space setting, the Particle MCMC algorithm used to sample from a generic target distribution, 
with static parameters = {^ji,^iiN) and latent state vector Xt = {xt,Ct,Vt,dt), {(^,xi-t\yi:t) proceeds 
by mimicking the marginal Metropolis-Hastings algorithm in which the acceptance probability is given by. 



Clearly, achieving this requires one to use a very particular structure for the proposal kernel in the MCMC 
algorithm. In particular the proposal kernel must take the form. 



Under this proposal, the acceptance probability for a standard Metropolis-Hastings algorithm takes the 
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The critical idea recognized in |14j in formulating the Particle MCMC algorithm is that the proposal 
distribution for vr {x'^.j\yi-f, can be sampled from via a Sequential Monte Carlo algorithm. 

6.1.1 Generic Particle MCMC Algorithm 

One iteration of the generic Particle MCMC algorithm proceeds as follows: 

1. Sample 6' r^q{e,-). 

2. Run an SMC algorithm with particles to obtain: 

N 



i=l 
t / N 

j=l \ i=l 



Wj \x. 



J 

Then sample a candidate path X[.^ ~ vf {xi-t\yi-t, ^')- 

3. Utilise the estimate of vf {yi:t\0) from the previous iteration of the Markov chain. 

4. Accept the proposed new Markov chain state comprised of {0\ X'^.^ with acceptance probability given 
by 

a {{6, X,,) , {e , X,,)) = mm ^1, ^ (^^^^|^) ^ ^ J (6-1) 

The key advantage of this approach is that the difficult problem of designing high dimensional proposals 
has been replaced with the simpler problem of designing low dimensional mutation kernels in the Sequential 
Monte Carlo algorithms embedded in the MCMC algorithm. In the paper fll|, this sampling approach in 
which Sequential Monte Carlo is used to approximate the marginal likelihood in the acceptance probability 
has been shown to have several theoretical convergence properties, see Section 4 of [14] for details. 

6.2 Adaptive Markov chain Monte Carlo (AdMCMC) Background 



There are several classes of adaptive MCMC algorithms, see [15]. The distinguishing feature of adaptive 
MCMC algorithms, compared to standard MCMC, is generation of the Markov chain via a sequence of 
transition kernels. Adaptive algorithms utilise a combination of time or state inhomogeneous proposal 
kernels. Each proposal in the sequence is allowed to depend on the past history of the Markov chain 
generated, resulting in many possible variants. 



Due to the inhomogeneity of the Markov kernel used in adaptive algorithms, it is particularly important 
to ensure the generated Markov chain is ergodic, with the appropriate stationary distribution. Several recent 
papers proposing theoretical conditions that must be satisfied to ensure ergodicity of adaptive algorithms 
include, [16], [E], [I7|, [E] and [19]. 

In [20] an adaptive Metropolis algorithm with proposal covariance adapted to the history of the Markov 
chain was developed. The original proof of ergodicity of the Markov chain under such an adaption was 
overly restrictive. It required a bounded state space and a uniformly ergodic Markov chain. In [15J a proof 
of ergodicity of adaptive MCMC under simpler conditions known as Diminishing Adaptation and Bounded 
Convergence is presented. In general it is non-trivial to develop adaption schemes which can be verified 
to satisfy these two conditions. In this paper we use the adaptive MCMC algorithm to learn the proposal 
distribution for the static parameters in our posterior In particular we work with an Adaptive Metropolis 
algorithm utilizing a mixture proposal kernel known to satisfy these two ergodicity conditions for unbounded 
state spaces and general classes of target posterior distribution, see [l5j for details. 



6.2.1 Adaptive Metropolis w^ithin Rao-Blackvi^ellised Particle MCMC 

In this section we present the specific details of the Adaptive Metropolis within Rao-Blackwellised Particle 
MCMC algorithm used to sample from the posterior on the path space of our latent factors and state 
space model parameters. This involves specifying the details of the proposal distribution in the PMCMC 
algorithm, 

q ((*, Xl:t, 0,..u Vl..t) , (6', x'l:U ^4'^)) = 9 (*', *) {x'l:t . ^V.t , Vi..t\Fl:t,T , C^..t,T,M,K ,^') • 

The proposal, q {^^^^~^\ for the static model parameters $ is given by an adaptive Metropolis proposal 
comprised of a mixture of Gaussians, one component of which has a covariance structure which is adaptively 
learnt on-line as the algorithm progressively explores the posterior distribution. The mixture proposal 
distribution for parameters $ is given at iteration j of the Markov chain by, 

q, •) = w,N i^; $(^-1), ^^S,^ + (1 - w,) N 1^; • (6-2) 

Here, Sj is the current empirical estimate of the covariance between the parameters of $ estimated using 
samples from the Particle Markov chain up to time j. The theoretical motivation for the choices of scale 
factors 2.38, 0.1 and dimension d are all provided in |15] and are based on optimality conditions presented 
in [21] and [22]. We note that the update of the covariance matrix, can be done recursively online via the 
following recursion, 

1 / , \ (6-3) 

Having specified the adaptive PMCMC proposal distribution for the static model parameters, we now focus 
on the proposal mechanism for the latent path space trajectories for the factors in the commodity model in 
the next subsection. 



6.3 Rao-Blackwellised Sequential Importance Sampling proposal 

The proposal kernel for the latent factors that will give us a Marginal Metropolis-Hastings framework is given 
by, vr (x'i:ti £,i-t^ (^'i f ^i:tWi:t,T-, ■ Sampling from this proposal can not be done exactly as it corresponds to 
sampling from the full conditional posterior distribution for the latent factors path genealogies. Instead we 
will utilise a specifically designed Sequential Monte Carlo algorithm to construct our proposal distribution. 

The SMC algorithm we develop here has a particular form which allows us to achieve two goals, the first 
that we can perform Rao-Blackwellisation via a Kalman-Filter for the latent factors Xi-.ti^i-.t- This optimal 
marginalised particle filtering approach and ideas are presented in [23] and similar ideas are discussed in 
|24| . Secondly, we develop a filter known as the Sequential Importance Sampling Resampling (SIR) filter for 
the latent factors Oi-t-, Vi-t- In doing this we can sample exactly from the prior for our latent factors (i.e. the 
Milstein discretization) and perform the importance weighting only dependent on the likelihood model. The 
algorithm for the Adaptive Metropolis within Rao-Blackwellised Particle MCMC is presented in Appendix 
[T2I Where the details of step 5. are also provided in Appendix [T2j The details of step 6. pertain to the 
Rao-Blackwellised SIR particle filter for constructing the proposal for the path spaces of the latent factors 
in the model. This requires some additional discussion and is presented below in Section 16. 31 

Remark 3: Development of a Rao-Blackwellized SIR particle filter proposal for the latent path trajec- 
tories proposal in the PMCMC algorithm, given the proposal static parameters $ has the advantage that 
it significantly reduces the variance in importance sampling weights. Achieved by exploiting the conditional 
linear Gaussian stat-space structure and a Kalman Filter. This will in turn increase the acceptance rate of 
our proposal mechanism. 

Remark 4: Utilising the SIR filter to construct a proposal for the latent factors 9i-t, Vi-t has an important 
advantage in the Milstein discretized state space model formulation we developed in that it can sample exactly 
from the prior for our latent factors and perform the importance weighting only dependent on the likelihood 
model. This is clearly beneficial in our setting since it allows us to avoid having to perform point-wise 
evaluation of the transition distributions given by /$ {Ot\Ot-i) and f^ {Vt\Vt-i) which can only be expressed 
as convolution integrals. 

The specific algorithm used to complete Step 7 of the Adaptive Metropolis within Rao-Blackwellised 
Particle MCMC algorithm is provided in Section [T2j Here we present the details for the sampling of a 
candidate proposal, conditional on the proposed state for the latent factor path spaces, xi-.t, Ci-.t, di-.t, Vi-.t 
as well as the estimation of the marginal likelihood p (-Fi:t,T|^') required for the calculation of the acceptance 
probability. This involves the development of a Rao-Blackwellised Kalman Filter scheme in conjunction with 
an SIR particle filter stage. 

The Rao-Blackwellisation via the Kalman Filter involves the following decomposition of the filtering 



distribution, 

TT {Xl:t,(.l:t, ^l.t, Vl;t\Fi:t,T , = VT (xi:t, Cl:t\Ol:t, Fl:t,T, *) 

X 7r{ei;t,Vi:t\F-i_:t,T,^) 

Next we note that the posterior distribution tt {xi:t,S,i:t\Oi:t,Vi:t, Fi:t,T,^) is exactly characterized by the 
Kalman Filtering estimate. We seek to exploit this fact in developing our sampling methodology. First we 
utilize a standard SIR particle filter with adaptive resampling to obtain the particle estimate, 

N 

yy 

i=l 

This particle estimate of the marginal then results in the joint posterior distribution estimate, 

N 



^{Xl:t,il:uOl..t,Vl..t\F^..t,T,^) = J^VF^^ (xi:t,ei:t|^S,^S,i^l:t,T,*) • 

i=l 



Now, clearly we can obtain vr (^i-t, S,i:t\Oi\, Vi^} , -fi:t,T, , for each particle i, an estimate of this distribu- 
tion, optimally via the Kalman Filter. In addition, it is easily shown that in this version of the SIR particle 
filter the unnormalised importance sampling weight for particle i is given by, see [23], 



We explain now the details of the Rao-Blackwellisation within the context of the SIR particle filter and 
then present the algorithm to perform this methodology. Given we have obtained a particle estimate 
l^i*!, , iy/*n,— 1- V at time t. We associate to each particle a Kalman Filter mean and covariance iif' = 

- ... „ ... ™ . .... k... ... _ . 

t\t~l ~ " n ^ " "n"n 

J-t -On-^t|f_iOn 

_J,T,,(«) I kT T/(«) 
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uTy{i) 



-'t\t ^t\t-l ^t\t-l"^ 

where = Gov (-F't,r|-P^i:t— i,t) and the predictive density is given by 

p iFt,T\Fl..t-l,T, Ol:t,Xl:t, ^l:t) = M (Ft,T] fl iF^..t-l,T^, O^-.t, Xl:t, ^l:t) , ^ 



Now that we have detailed each aspect of the model and each aspect of the proposed sampling algorithm 
to allow us to calibrate and estimate the Bayesian commodity model developed, we will now present the 
results and analysis. This will be done in two stages, first we present details of a comprehensive synthetic 
simulation studv, followed bv a real data analvsis. 



7 Results and Discussion 



In this section we present studies on the performance of the Bayesian model formulation, the Milstein 
discretization framework and the Adaptive PMCMC sampling framework under Rao-Blackwellization. To 
achieve this we first consider a case study consisting of synthetic data analysis. After demonstrating the ac- 
curacy of the estimation and calibration methodology we developed, we perform a detailed analysis utilising 
actual data representing a panel of crude oil futures. 

7.1 Synthetic Data Example 

In this section we present detailed analysis of the performance of the Adaptive PMCMC algorithm in the 
context of the Bayesian state space model non-linear filtering and parameter estimation framework that 
we have developed. This study aims to systematically assess the following model and AdPMCMC sampler 
aspects: the ability to estimate static model parameters (perform calibration) in high observation noise; the 
impact of parameter estimation under vague priors for static parameters corresponding to drift and volatility 
of each of the four discretized s.d.e. factor models; the impact of increasing the number of futures contracts 
with different maturities on parameter and latent factor estimation; and finally, the impact of the Milstein 
scheme truncation for p. 

7.1.1 Analysis of Estimation Accuracy versus Number of Futures Contracts 

The simulation study for this section proceeded by first sampling a set of static model parameters from 
the priors which we will utilise to generate synthetic data. We will then obtain posterior estimates of these 
parameters given the generated data. The 'true' static parameters utilised in this study, were given by = 
{f3, fi(^, K^, fiy, Kv,cr^,a^,av,^F) = ^rn = (0.2,0.1,0.4,0.2,0.2,0.5,0.5,0.5,4) and seasonality components 
wiai = 0. These parameter settings were then utilized to generate a set of latent process trajectories which 
were then in turn utilized to generate a set of futures curves with observations in noise. This produced a set of 
synthetic latent paths, for 100 days T = 100, for xi-.T, Ci-.T, di-.T, Vi-.T which are presented in the upper panel 
of Figure El The lower panel of Figure H] presents the resulting log spot price each day. We then utilised these 
simulated trajectories to generate a panel of futures curves according to the observation equations developed 
in Section We utilised a panel of futures curves corresponding to 10 futures contracts, each separated by 
maturities of 30 days, the first contract rolling over in 29 days time, followed by the next contract maturing 
in 59 days, all the way up to the tenth contract which at time t = 1 of the observations still had 299 
days remaining until maturity. Using these latent process, the noisy observations of the futures curve was 
generated over 100 days of observations using T,f = diag (^api^j,^y . . . ,(t|,^^^^^ = (1, . . . , 1). Additionally, 
unless otherwise stated we ran the AdPMCMC sampler presented in Algorithm 1 for J = 200, 000 Markov 
chain samples with = 500 particles. The discard Markov chain samples for burn-in was selected as 
50, 000. The discretization interval considered was daily, At = 1 and the p-th order truncation of the 
Milstein discretization used to calculate Or> involved setting; » = 100. 



In this section, we consider uniform priors on the static state space model parameters, corresponding to 
hyper-parameter settings specified according to ~ f7 (0, 10), cj^ ~ U{0, 10), /x^ ~ (0, 10), ^ U (0, 10) 
and cj^ ~ C/(0, 10). Additionally, we assume the covariance matrix of the additive multivariate Gaussian 
observation noise on the log of the futures contracts over different maturities, on any given day t, is diagonal. 
We then consider the observation noise variance for the Tj-th futures contract on day t, denoted cr^. utilised 
to generate the observations as given by S^r = 41 which is corresponding to a coefficient of variation in the 
noise of around 20%. Note, in the Bayesian model we treat this observation noise variances as unknown 
quantities a priori that must be estimated jointly with the other model parameters from the posterior. 

We consider a panel of N different futures contracts with maturities Ti , . . . , T/v and futures prices at 
time t denoted by -Ft,Ti, • • • ,Ft^TM- Note, the filtering framework is set up such that the latent process is 
the real process with real parameters used and the observation process is the real process with risk neutral 
parameters used. In the following section we analyse the ability to perform jointly parameter estimation 
(calibration) and filtering for the latent states of the model as we vary the number of contracts observed on 
the futures curve on any given day from G {1, 5, 10}. 

In the following simulation results we demonstrate the accuracy of the joint calibration (parameter 
estimation) and filtering (state estimation) using the posterior distribution derived in Equation 15.31 and the 
Rao-Blackwellised Adaptive Particle MCMC algorithm. The first results presented are for the parameter 
estimations as a function of the number of Futures contracts considered, where we consider N £ {1, 5, 10}. 
In Figure [1] we present the MMSE and 95% posterior confidence intervals for each model parameters real 
and risk neutral. Each plot contains the parameter estimates as well as the true parameter values used 
to generate the data represented by circles. The top panel contains the posterior estimates obtained for 
a panel of 100 days of data with = 1 futures contract, the middle panel has posterior estimates for a 
panel of 100 days of data with = 5 futures contracts and the lower panel has these results for 100 days 
of observations with A^ = 10 contracts. Each contract considered is separated by 30 days so we have 30 
day contract, 60 day, through to maturity TiO which is a 300 day contract. The results demonstrate that 
clearly the parameter estimation accuracy is improving as we add more contracts to the panel of data, with 
A^ = 10 contracts having very accurate parameter estimates. 

In Figure [2] we present the trace plots of the model parameters from the Rao-Blackwellised Adaptive 
PMCMC sampler for futures contract panel data contains A^ = 10 contracts separated by 30 day maturities 
over 100 trading days. We have separated the results per factor and present the real world and risk neutral 
parameters on the same trace plot. These results are for uninformative priors and we see reasonable mixing 
for the Rao-Blackwellised AdPMCMC sampler updating in dimensions 4T + 18 each iteration. 

In Figure [3] we present the Kalman filtered and particle filtered estimates of the latent path space for each 
of the four factors for short and long term dynamics and the two stochastic volatility terms, for the case 
of A^ = 10 contracts separated by 30 day maturities over 100 trading days in the panel of futures curves. 
The results demonstrate that the short term factor and volatility terms are very accurately estimated. The 
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Figure 1: Filled circles are the estimated parameter MMSE and error bars represent 95% posterior confidence 
intervals. Open circles are true parameter values used to generate the data. Parameters are presented in 
order of real world (left) and risk neutral (right). 
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Figure 2: Trace plots of the simulated Markov chain state trajectories generated by the Rao-Blackwellised 
AdPMCMC sampler for the model parameters of the multi-factor s.d.e. model. 

long term factor estimates the mean reversion level well, and includes the true latent process within a 95% 
posterior confidence interval. 
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Figure 3: True latent factor model s.d.e. trajectories (solid line) versus Rao-Blackwellised AdPMCMC 
sampler filter estimates for the MMSE (dashed line) and 95% posterior confidence intervals (dotted lines) 
for each factor: short term and long term trends and volatility. The Milstein SIR proposal utilised p = 100 
and 500 particles. 

Finally, the most important result in this analysis is a comparison of the estimated log spot price versus the 
estimated MMSE log spot price and 95% posterior confidence interval, see Figure HI As above, this example 
corresponds to the situation in which the futures contract panel data contains 10 contracts separated by 
30 day maturities over 100 trading days of observations. The Milstein SIR proposal considered in the 
Rao-Blackwellized Adaptive PMCMC algorithm utilised truncation of p = 100 and 500 particles. 

We also note that simulation studies were also performed utilising futures contracts generated with 
seasonal components, the resulting estimation accuracy were equivalent to those presented. Additionally, 
we varied the signal to noise ratio by decreasing observation error to ap = 1 and observed a corresponding 
increase in accuracy of the estimation. 

7.1.2 Actual Data - Oil Commodity Futures 

In this section we take a panel 10 futures contracts with daily price data corresponding to the contracts for 
crude light oil commodity futures over 150 trading days from 2/01/1998 to 6/08/1998. The maturities of 
the 10 contracts are 21 days for the first contracts initiation and each contract there after has a maturity 
increasing by 20 days. Therefore the last contract has a maturity from its initiation of 200 days. We 
estimated the posterior model MMSE parameters and posterior confidence intervals for the developed in 
Section [5] using this panel of oil futures data. This was achieved via the Rao-Blackwellised AdPMCMC 
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Figure 4: True latent log spot price trajectory versus Rao-Blackwellised AdPMCMC sampler filter estimated 
log spot price MMSE and 95% posterior confidence interval. 

algorithm for the four factor model developed, again with a calibration and filtering formulation and the 
100,000 iterations of the sampler. The average acceptance rate of the PMCMC sampler developed was 
28% and 500 particles were utilised. The Milstein SIR proposal utilised p = 100 and the calibrated model 
parameters obtained are presented in Figure [5] along with 3 posterior standard deviations of these estimated 
parameters using the model developed in Section El 




Figure 5: Calibrated MMSE posterior model parameters and posterior confidence intervals. This example 
corresponds to the situation in which the crude oil futures contract panel data contains 10 contracts separated 
by 20 day maturities over 150 trading days. The Milstein SIR proposal utilised p = 100 and 500 particles. 

In Figure [6] we present the latent path space trajectory estimates for the four factors over 150 days fitted 
to the panel of oil futures data. The plots contain the posterior MMSE estimates for each factor and the 
95% posterior C.I. estimates. In addition, we present the estiamted log spot price for the crude oil futures 
contract data, based on 10 contracts with 150 days of data in the panel. 

Next we present the MMSEP and the associated posterior predictive 95% confidence intervals obtained 
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Figure 6: Panel 1: Short term dynamic Xi-.T, Panel 2: Long term dynamic Ci:Ti Panel 3: Volatility term 
6i;T, Panel 4: Volatility Vi:T, Panel 5: Log spot price estimated Si:T- Solid line is filtered MMSE and gray 
band is posterior 95% C.L 

from the predictive distribution of the log futures curves at some time point r defined by Equation 17.11 

P{Pt,Ti , Pt,T2^ • ■ • : Pt.Tn \Pi:T,Ti , Pl:T,T2 7 • • • 7 ^1:T,Tn) 
N 

= Yl d^i^-^.TjXrXr, Or)TT (*, Xl:T, Cl:T, ^'l:T, Vi;t\Fi;T,Ti , ^'l:T,T2 , • ■ • , Fl.,T,Ti^) d^dxi:T.TdCl:T,Td01:T,TdVl;T.T 
•' n=l 

(7.1) 

Results in Figure [7] represent from top to bottom in the figure panels on the left the within sample and 
on the rieht the out of samnle forecast loe futures nrices after integrating out uncertainty in the latent 



process states and parameters, resulting from the posterior predictive distribution. The shaded gray area 
corresponds to the predicted 95% posterior predictive confidence interval for the log of the futures prices for 
contracts 1, 5, 10 observed over 150 days and predicted out of sample over an additional 5 days. The first 
contract in the top row of results is the contract closest to maturing, the 5th contract in the middle row of 
results is contracts maturing at 100 days and the 10-th contract in the bottom row of results are those with 
maturties of 200 days. The solid line represents the observed time series of daily log futures price at close 
of market over 150 days within sample and 5 days out of sample. The dashed line represents the estimated 
posterior MMSE and predicted MMSE of the log futures price. Clearly we see that the observed futures 
price curves are contained within the 95% posterior C.I of the estimated model at all times, and secondly 
that the estimated MMSE and forecast MMSE out of sample indicate acturate calibration and estimation 
of the commodity model on the real crude oil futures data. 




Figure 7: Left Panels: within sample prediction, Right Panels: out of sample forecasts. Top Row: contracts 
with 20 days maturity. Middle Row: contracts with 100 days maturity. Bottom Row: contracts with 200 
days maturity. This example corresponds to the situation in which the crude oil futures contract panel data 
contains 10 contracts separated by 20 day maturities over 150 trading days. Solid line observed log futures 
price, dashed line predicted MMSE log futures price, gray shading 95% posterior C.I. 



8 Conclusions 



In this paper we have developed a novel general multi-factor model for commodity spot prices and futures 
valuation. In particular we extend the multi-factor model long-short stochastic differential equation (s.d.e.) 



model developed in [1] and [2] in two important aspects: firstly we allow for both the long and short term 
dynamic factors to be mean reverting and secondly we develop an additive structural seasonality model. We 
developed a non-trivial state space model representation for such a model and in doing so we also extend 
the literature in this area by allowing for futures and options contracts in the observation equation of the 
Milstein discretized non-linear stochastic volatility state space model. 

Next to perform estimation and calibration under this model we derived a closed form solution for the 
futures prices and then developed a novel numerical methodology based on Sequential Monte Carlo to 
perform calibration of the model jointly with the filtering of the latent processes for the long-short and 
volatility dynamics. In this regard we explore and develop a novel methodology based on an adaptive Rao- 
Blackwellised version of the Particle Markov chain Monte Carlo methodology for the joint calibration and 
filtering. In doing this we deal accurately with the non-linearities introduced into the filtering framework. 
We perform analysis on real data for oil commodities. 
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9 Appendix 1: Proof of Futures Price 



Proof: Here we present the proof of Theorem 1 for the futures price from the real latent process s.d.e. 
model. To proceed we obtain the risk neutral process from the real process we need to deduct the risk 
free premium from the long term trend to account for the market price of risk where the risk premium is 
denoted by X{x,t). Next, we denote the Futures contract price at time t with maturity of the contract at 
ThyFit,T). 

F{t, T) = E* [exp (X(r)) \X{t)] = E* [exp (/(T) + + Xt + ^t) 1/(0, Ct*, Xt* , ] , 

where E* denotes the expectation taken with respect to the risk neutral process. To solve for this ex- 
pectation, we must first derive the expression for the transition density of the state variables, denoted 
by P ('^T + Xt + ^T' -^l^t' Xt 1 ) — P i^TlYtjt). This is obtained for the continuous s.d.e. latent 
process model by writing down the either the Kolmogorov-backward or forward equation. Here we work 
with the backward equation (dropping the * notation on state variables for convenience) with the condition 
that, p {YT,T\Yt,t = T) = 5{Yt — Yf). Next, to obtain the futures price we multiply each term in the 
Kolmogorov-backward equation by exp (Xf) and then integrate with respect to Xt- Note, the notation Ff^T 
refers to F(Yt,t = T,T). This allows us to express each term in the pdc with respect to the futures price. 



t,T 



^dxtdit dxtdVt 
with condition F{Yt, t = T,T) = exp(XT). 

Next, since the state variables enter into the above pde in a linear fashion, therefore the solution belongs 
to the class of exponentially affine models. Hence, we can look for a solution of the form, 

Ft,T = exp (So(r) + 5i(r)6 + 52(r)xt + B3{T)et + B^{T)Vt) . 

Hence, substitution of this solution results in a system of 5 ODE's which after substitution and grouping of 
like terms, then divide both sides by Ft^T and group terms according to latent factors to obtain: 

dB2{T) 



dT 
dT 



-l3*B2iT); 
0; 



= -X^Bsir) - K*yB4T) + \a\ {B^irf + pevavB^(T)B^{T) + \bI{t); 

UT Z Z 

= -^^1 [Bi{r)f - \al [B,{r)f - B,{r)f,l + X,B,{r) - „*yB,{r) - p^^a^a^B,{T)B,{r) 

Note, the condition that F{Yt,t = T,T) = exp(Xj) needs to be satisfied, in addition we assumed the form 
of the risk premium X{Xt,t) = A4\A4 = — ^VVt (this allows us to remove the volatility Vt from the futures 
orice. however still keeo it in the ootions orice). 



The solution to this system of ODE's is obtained by first applying the constraint at time T given by 
F{Yt,t = T,T) = exp(XT) which automatically sets Bi{0) = ^2(0) = ^3(0) = 1 and Bo{0) = ^4(0) = 0. 
We can then obtain solutions B2{t) = exp(— /3*r), Bi{t) = exp(— K|r) and with initial condition ^3(0) = 1 
we get Bs{t) = 1. Finally, substituting the solution for B^It) into the expression for '^^^^^ and setting 
A4 = — ^ results in ODE 

^^M. = [p^yav - K*y] B^{t) + \al \B^{T)f 

with the initial condition that -84(0) = 0. This is recognized as Ricatti's equation of the 2nd kind one can 
obtain a solution to a transformed linear 2nd order ode and then in turn obtain a solution to the Ricatti 
equation which produces a solution .64 (r) = 0. Finally, substitution of these solutions for Bi{t),B2{t) and 
B4{t) into the expression for results in an ode which can be integrated (w.r.t. r) explicitly to obtain 

a solution for Bo{t) as follows: 

Mr) = -\4f exp(-2K|r)dr - j exp(-2rr)dr- 

j exp (— K|r)/x|cir + Ao j exp{—P*T)dT — p^-^a^a^ J exp(— ^*t) exp(— K|r)(ir 

= ^ exp(-2K|r) + ^ exp(-2/3V)+ 

^ exp(-K^r) - — exp(-^ r) + ^^r^ exp(-^ r - k^t) 

□ 

10 Appendix 2: Proof of Strong Stochastic Taylor Representation of 
s.d.e. Factors 

The strong Taylor stochastic expansion known as the Milstein scheme in a multivariate setting is presented 
below. We present the i-th component of a general n-dimensional s.d.e., with m-dimensional Weiner process 
as, 

m 

dXl = 0} (t, Xt) dt + Y^ U'^ (t, Xt) dWi- 

Throughout the remainder of this section wc will adopt here the operator notation of Klocdcn and Platen 
(1999). In their book, p. 346 they demonstrate that the i-th component of the Milstein scheme under the 
Ito integrals is then given by, 

m 

Xl = X\_^ ^a\t- 1, At + ^ 6*'^' [t - 1, Xt-x) l\Wl_^ 

m 

where, LP = J2^=i ^^'"'^ Ito multiple integral I(j^^j^)At is given by, 

hn,n)M= / / dWlldWil. 

*i tn. *^ tn. 



These integrals have the useful properties that, 

and I(j-i^j^)At with ji 7^ j2 not easily expressed. 

As pointed out in Kloeden and Platen (1999), when ji ^ j2 the Ito and Stratonovich integrals are equal, 



It will be easier to obtain the approximation under the Stratonovich integrals which under a p-th order 
truncation is given by, 

'^{ji,j2)At ~ A^ (^2'^jiCj2 + ^/Pp {fJ'ji,pCj2 ~ Mj2,p0l)^ 



r=l 



where Q, pj^p, vj^r and tpj^r are all independent AA(0; 1) Gaussian random variables with, 



= 1- J_v-1 

PP 12 27r2^r2' 

r=l 



1 



We note that a bivariate approximation for the Heston model was studied in Stump (...) and a rec- 
ommendation for p is also provided in Kloeden and Platen (1999), who suggest p > for some positive 
constant K. 

Hence, the bivariate s.d.e.'s for 9t and Vt given by, 



dVt = {nv — Vi) cii + av^/VtdZv 



will first be recast with respect to independent Weiner processes dW\ and dW2 as follows 

dOt = \/VtdWi 

dVt = (pv - KyVt) dt + ayy/Vt ^pvedWi + y^l - ^^^^14^2^ • 



This will result in the following specifications, (n=m=2): 

{t, et, Vt) = 0; {t, Ot, Vt) = {fiy - i^vVt) ; b^'^ {t, Ot, Vt) = h"^^ {t, Ot, Vt) = 
fe^'i {t, 9t, Vt) = avVVtPve; {t, et, Vt) = avVVt^l - P^g 

L^h"^^ it, Bt, Vt) = h"'^ {t, 9t, Vt) ^6^'^ (t, Ot, Vt) + b^'^ (t, 9t, V) ^6^'^ (t, Ot, V) = -aypve 

lH^'' {t, et, Vt) = b^'^ it, et, Vt) -^b"'^ it, et, Vt) + 1^-^ it, Ot, Vt) -^b''' it, et, Vt) = \ov^Jl^e 

L'b^'^ it, Bt, Vt) = 0; LH^'^ it, Ot, Vt) = 0; 

L^t^'l it, et, Vt) = b''^ it, Bt, Vt) ^62,1 y^^ + ^2,2 y^^ ^^2,1 y^ ^ ^^lpy^^i_p2^^ 

L'b^^^ it, Bt, Vt) = b''' it, Bt, Vt) ^6^,2 (^^ y^^ + ^2,1 (^^ y^ ^^2,2 (^^ y^^ ^ ^^p^., ^1 - p^., 

lH^'^ it, Bt, Vt) = b^'^ it, Bt, Vt) it, Bt, Vt) + 62.2 y^^ _d_^2,2 y^-^ ^ 1 ^2 _ ^2^^) 

L^b^'i = 6i'i it, Bt, Vt) ^62,1 y^) + 52,1 ^^2,1 ^ 1 ^2^^2^^ 

This leads to the following bivariate strong Milstein discretization scheme, 

et = et-i + ^/Vt-iAtne^t + -2'^vPveI{i,i)At + ~ Py6t-f(2,l)At 

Vt = Vt-i + ipv - i^vVt-\) At + avPveV^t-iAtng^t + crv\JVt-i (l - Pyg) Atnv,t 

]_ 2 2 ^2 / 2 ^ 2 / 2 ^ 2 / 2 \ 

+ 2^vPye-^(i,i)At + -^'^vPvey/'^- Pye-^(i,2)At + -^'^vPvey'^ - Pve^{2,i)At + 2<^v (1 - Pye) ^(2,2)At 

« Vt-i + (^v- - KyV^-i) At + fTv/CyeVVt-iAtne^f + (Tv\JVt-i (l - /Cyg) Atny^t 

+ \(^vPve i^tn^t - Ai) + ^c^YPve^i - PveJ(i,2)At + ^^^yP^^e-v/l - Pye-^(2,i)At + (l - Pve) {^tnl^t - 

where ng^t <^iid ny^t are i.i.d. standard normal random variables and we will approximate the mixed ex- 
pressions J2,iiB) and Ji,2(^) with p-th order truncations of the series, when we perform simulation of the 
series. 

11 Appendix 3: Simulating a Strong Stochastic Taylor Scheme 



Simulation of the strong stochastic Taylor scheme for a Miltstein approximation is give next. This will form 
an important aspect of the proposal mechanism in the SIR fitler contained in the Particle MCMC algorithm. 
We wish to simulate new realizations for Bt and Vt given Bt-i and Vt-i and model parameters at time t, 



Simulation of a bivariate Milstein Scheme. 



1. Initialisation: 



i. Specify time discretization At and truncation p. 



2. Simulation: 



Fori = l,...,N 

i. simulate ny\_-^ ' 

ii. simulate riy ~ 
to obtain {riy^ s 

iii. simulate Ci\C. 




For r = 1, . . . ,p 

iv. simulate -01*]. > '02*r > ^i*r' ^2*r- each independently from M (0, f ). 
3. Evaluate: 

For i = l,...,N 

i. evaluate ■^i2i^) using equation [10] and simulated terms from stage 2. 

ii. evaluate Milstein approximations to obtain new samples 9^^^ and V^^^ given O^^^ and V^^-^ 

Table 1: Generation from the 2-d Milstein strong stochastic Taylor scheme, truncated at p-th order. 



12 Appendix 4: Adaptive Metropolis within Rao-Blackwellised Particle 
MCMC Algorithm 



Step 2: Sample the static model parameters. 



if ui > wi then 



/* perform an adaptive random walk move */ 




else 



/* perform a non-adaptive random walk move */ 




Algorithm 1: Adaptive Metropolis within Rao-Blackwellised Particle MCMC 



Input: Initial Markov chain state , xS > > > ^i?) • 

Output: J Markov chain samples {*(^\ xS' ^S, ^b?}j=i:J ~ P (*, Xi:t, 6:t, ^i:t, ^i:t|-Fi:t,T). 
begin 

1. Set initial state x^j , , , ^i-? j deterministically or by sampling the priors, 

repeat 

2. Sample given from an adaptive multivariate Gaussian Metropolis proposal, to obtain 
proposal 

3. Sample Xi-.t, ^i-.t, ^i-.t, Vi-.t via a combination of Rao-Blackwellizing Kalman Filter and SIR particle 
filtering, conditional on to obtain proposal {^',Xi:t^^i:ti^i:ti^(:t) 

4. Accept proposal ,Xi:t'^i:t^^i:t'^i:t) with Metropolis-Hastings acceptance probability given by 

« ((*^^-^ ^!fr'\ ^r'O ' ' ^^*)) 

= min i ^ 

and increment j = j + 1. 



until j=J 

end 



Algorithm 2: Step 3: Particle Filter Proposal - Conditional on the sampled static model parameters 
6' , run an SIR particle filter with N particles 

3(i). Initialise the SIR particle filter for Xq = {Oq^Vq) with N particles {^6'^*^ V^o'^)}i=i:Ar 

3(ii). Initialise the Kalman filter mean and covariance for Xo^^o*'' denoted by yit^*'' = ^//^*''(0),//^*''(0)j and 

^* ~^(x,O(0) 
repeat 

3(iii). Perform Kalman filter update for the N particles at iteration t — 1, to update /i^!^^ and S^^^^ 
to ixf' and E^*^ according to recursions in (Section 6.3) 

3(iv). Perform Mutation of the N particles at iteration t — 1 to obtain new particles at iteration t via 
prior, Xj'^ ~ M [xt; ft {xt-i,0') , af). Append the new state at time t to the previous particle 
geneology to obtain Xq.^ ion i = 1, . . . , N 

3(v). Evaluate the importance sampling weight for the A'^ particles as, z-th weight is 
until t=T ; 

3(vi). Normalise the importance sampling wei ghts Wi ' = ;^'~(,) 

3(vii). If the Effective Sample size is less than 80% then resample the particles at time t using stratified 
resampling based on the empirical distribution constructed from the importance weights to obtain 
new particles with equal weight. 

3(viii). Sample one proposed path trajectory from the particle estimated proposal distribution given by 

X[.j, ~ P (rri:T|yi:T, 9') = Eili W^lS (xi:t) 

•''liT 

3(ix). Evaluate the marginal likelihood given by p{yi:TW)- 



